

% function Print_Equ_NK

%% Declare Variables, Parameters, and Steady States
% Variables
syms        Eta_Z Eta_Zp Eta_Z_N Eta_Z_Np Eta_Z_H Eta_Z_Hp ...
            Eta_M_dom Eta_M_domp Eta_M_ext Eta_M_extp ...
            Eta_Y_H Eta_Y_Hp Eta_p_F Eta_p_Fp ...
            N_N N_H w_N w_H mc_N mc_H Y_N Y_H Profit_N Profit_H PAC_N PAC_H ...
            Y_Np Y_Hp ...
            L_N L_H C_N C_T C_H C_F ...
            p_T p_N p_H p_F Pir Pirp Pir_N Pir_H Pir_F ...
            Lag_p_T Lag_p_N Lag_p_H Lag_p_F Lag_p_Tp Lag_p_Np Lag_p_Hp Lag_p_Fp ...
            p_Np p_Hp Pir_Np Pir_Hp ...
            Bp_dom ir_dom ir_ext rr_dom rr_ext ...
            ir TAU B T G ir_s Pir_H_s Pir_F_s dE dEp Y dY C Profit PAC ...
            Lag_B Lag_Bp Lag_ir Lag_irp Lag_Eta_p_F Lag_Eta_p_Fp Lag_Y Lag_Yp;
varlist=   [Eta_Z Eta_Zp Eta_Z_N Eta_Z_Np Eta_Z_H Eta_Z_Hp ...
            Eta_M_dom Eta_M_domp Eta_M_ext Eta_M_extp ...
            Eta_Y_H Eta_Y_Hp Eta_p_F Eta_p_Fp ...
            N_N N_H w_N w_H mc_N mc_H Y_N Y_H Profit_N Profit_H PAC_N PAC_H ...
            Y_Np Y_Hp ...
            L_N L_H C_N C_T C_H C_F ...
            p_T p_N p_H p_F Pir Pirp Pir_N Pir_H Pir_F ...
            Lag_p_T Lag_p_N Lag_p_H Lag_p_F Lag_p_Tp Lag_p_Np Lag_p_Hp Lag_p_Fp ...
            p_Np p_Hp Pir_Np Pir_Hp ...
            Bp_dom ir_dom ir_ext rr_dom rr_ext ...
            ir TAU B T G ir_s Pir_H_s Pir_F_s dE dEp Y dY C Profit PAC ...
            Lag_B Lag_Bp Lag_ir Lag_irp Lag_Eta_p_F Lag_Eta_p_Fp Lag_Y Lag_Yp];
% Parameters
syms        PP_RHO_Z PP_RHO_Z_N PP_RHO_Z_H ...
            PP_RHO_M_dom PP_RHO_M_ext PP_RHO_Y_H PP_RHO_p_F ...
            PP_ALPHA_N PP_ALPHA_H PP_EPSILON_N PP_EPSILON_H ...
            PP_THETA_N PP_THETA_H ...
            PP_Cons_FracT PP_Cons_FracH PP_Cons_ElasTN PP_Cons_ElasHF ...
            PP_Taylor_Pi PP_Taylor_Pi_H PP_Taylor_Pi_N PP_Taylor_dE PP_Taylor_dY PP_Taylor_ir...
            PP_Flag_dE PP_Flag_Price ...
            PP_FiscalFlag_TAU PP_FiscalFlag_T PP_FiscalFlag_G PP_FiscalFlag_B PP_FiscalFlag_B_Real;
pplist=    [PP_RHO_Z PP_RHO_Z_N PP_RHO_Z_H ...
            PP_RHO_M_dom PP_RHO_M_ext PP_RHO_Y_H PP_RHO_p_F ...
            PP_ALPHA_N PP_ALPHA_H PP_EPSILON_N PP_EPSILON_H ...
            PP_THETA_N PP_THETA_H ...
            PP_Cons_FracT PP_Cons_FracH PP_Cons_ElasTN PP_Cons_ElasHF ...
            PP_Taylor_Pi PP_Taylor_Pi_H PP_Taylor_Pi_N PP_Taylor_dE PP_Taylor_dY PP_Taylor_ir ...
            PP_Flag_dE PP_Flag_Price ...
            PP_FiscalFlag_TAU PP_FiscalFlag_T PP_FiscalFlag_G PP_FiscalFlag_B PP_FiscalFlag_B_Real];
% Steady States
syms        SS_Z_N SS_Z_H ...
            SS_N_N SS_N_H SS_w_N SS_w_H SS_mc_N SS_mc_H SS_Y_N SS_Y_H ...
            SS_Profit_N SS_Profit_H SS_PAC_N SS_PAC_H ...
            SS_C_T SS_C_N SS_C_H SS_C_F SS_L_N SS_L_H ...
            SS_p_T SS_p_N SS_p_H SS_p_F SS_Pir SS_Pir_N SS_Pir_H SS_Pir_F ...
            SS_Bp_dom SS_ir_dom SS_ir_ext ...
            SS_ir SS_TAU SS_B SS_T SS_Y SS_C SS_Profit SS_PAC SS_G SS_C_H_s;
sslist=    [SS_Z_N SS_Z_H ...
            SS_N_N SS_N_H SS_w_N SS_w_H SS_mc_N SS_mc_H SS_Y_N SS_Y_H ...
            SS_Profit_N SS_Profit_H SS_PAC_N SS_PAC_H ...
            SS_C_T SS_C_N SS_C_H SS_C_F SS_L_N SS_L_H ...
            SS_p_T SS_p_N SS_p_H SS_p_F SS_Pir SS_Pir_N SS_Pir_H SS_Pir_F ...
            SS_Bp_dom SS_ir_dom SS_ir_ext ...
            SS_ir SS_TAU SS_B SS_T SS_Y SS_C SS_Profit SS_PAC SS_G  SS_C_H_s];
        

%% Equations
%--------------------------------------------------------------------------
% Note: the aggregate quantities Y,C,G,T,Bh,Bg,PAC,Profit
%       are all normalized by SS_w*SS_N. w and N are in their original
%       levels.
%--------------------------------------------------------------------------
e       =   sym('e%d',[56,1]);
%--------------------------------------------------------------------------
% Dynamics of Exogenous Variable (6)
%--------------------------------------------------------------------------
id      =   1;
e(id)   =   -Eta_Zp+PP_RHO_Z*Eta_Z;
id      =   id+1;
e(id)   =   -Eta_Z_Np+PP_RHO_Z_N*Eta_Z_N;
id      =   id+1;
e(id)   =   -Eta_Z_Hp+PP_RHO_Z_H*Eta_Z_H;
id      =   id+1;
e(id)   =   -Eta_M_domp+PP_RHO_M_dom*Eta_M_dom;
id      =   id+1;
e(id)   =   -Eta_M_extp+PP_RHO_M_ext*Eta_M_ext;
id      =   id+1;
e(id)   =   -Eta_Y_Hp+PP_RHO_Y_H*Eta_Y_H;
id      =   id+1;
e(id)   =   -Eta_p_Fp+PP_RHO_p_F*Eta_p_F;
%--------------------------------------------------------------------------
% Firms' Decisions (10)
%--------------------------------------------------------------------------
% Production Function (2)
id      =   id+1;
e(id)   =   -SS_Y_N*exp(Y_N) ...
            + SS_Z_N*exp(Eta_Z+Eta_Z_N)*( SS_N_N*exp(N_N) )^(1-PP_ALPHA_N);
id      =   id+1;
e(id)   =   -SS_Y_H*exp(Y_H) ...
            + SS_Z_H*exp(Eta_Z+Eta_Z_H)*( SS_N_H*exp(N_H) )^(1-PP_ALPHA_H);
% Marginal Cost (2)
id      =   id+1;
e(id)   =   -SS_mc_N*exp(mc_N) ...
            +(SS_N_N*SS_w_N)*exp(N_N+w_N)/(SS_Y_N*exp(Y_N))/(1-PP_ALPHA_N);
id      =   id+1;
e(id)   =   -SS_mc_H*exp(mc_H) ...
            +(SS_N_H*SS_w_H)*exp(N_H+w_H)/(SS_Y_H*exp(Y_H))/(1-PP_ALPHA_H);       
% Price Setting (Phillips Curves) (2)
id      =   id+1;
e(id)   =   -(Pir_N ) ...
            +1/(1+SS_Pir_N)*PP_EPSILON_N/PP_THETA_N*( mc_N-p_N ) ...
            +(1+SS_Pir_N)/(1+SS_ir)*(Pir_Np);
id      =   id+1;
e(id)   =   -(Pir_H - PP_Flag_Price*(1+SS_Pir_H)*dE) ...
            +1/(1+SS_Pir_H)*PP_EPSILON_H/PP_THETA_H*( mc_H-p_H ) ...
            +(1+SS_Pir_H)/(1+SS_ir)*( Pir_Hp - PP_Flag_Price*(1+SS_Pir_H)*dEp );
% Price Adjustment Costs (2)
id      =   id+1;
e(id)   =   -( PAC_N ) ...
            + 0;
id      =   id+1;
e(id)   =   -( PAC_H ) ...
            + 0;
% Profit (2)
id      =   id+1;
e(id)   =   -( SS_Profit_N+Profit_N ) ...
            + SS_p_N*SS_Y_N*exp(p_N+Y_N)-exp(N_N+w_N)*SS_N_N*SS_w_N- ( SS_PAC_N+PAC_N );
id      =   id+1;
e(id)   =   -( SS_Profit_H+Profit_H ) ...
            + SS_p_H*SS_Y_H*exp(p_H+Y_H)-exp(N_H+w_H)*SS_N_H*SS_w_H- ( SS_PAC_H+PAC_H );

%--------------------------------------------------------------------------
% Households' Decisions (9)
%--------------------------------------------------------------------------
% Consumption Basket (4)
id      =   id+1;
e(id)   =   -SS_C_T*exp(C_T) ...
            +SS_C*exp(C)*PP_Cons_FracT*( SS_p_T*exp(p_T) )^(-PP_Cons_ElasTN);
id      =   id+1;
e(id)   =   -SS_C_N*exp(C_N) ...
            +SS_C*exp(C)*(1-PP_Cons_FracT)*( SS_p_N*exp(p_N) )^(-PP_Cons_ElasTN);
id      =   id+1;
e(id)   =   -SS_C_H*exp(C_H) ...
            +SS_C_T*exp(C_T)*PP_Cons_FracH*( SS_p_H/SS_p_T*exp(p_H-p_T) )^(-PP_Cons_ElasHF);
id      =   id+1;
e(id)   =   -SS_C_F*exp(C_F) ...
            +SS_C_T*exp(C_T)*(1-PP_Cons_FracH)*( SS_p_F/SS_p_T*exp(p_F-p_T) )^(-PP_Cons_ElasHF);
% Price Index (2)
id      =   id+1;
% e(id)   =   -1 ...
%             +PP_Cons_FracT*( SS_p_T*exp(p_T) )^(1-PP_Cons_ElasTN) ...
%             +(1-PP_Cons_FracT)*( SS_p_N*exp(p_N) )^(1-PP_Cons_ElasTN);
e(id)   =   -0 ...
            +PP_Cons_FracT*( p_T )+(1-PP_Cons_FracT)*( p_N );
id      =   id+1;
% e(id)   =   -1 ...
%             +PP_Cons_FracH*( SS_p_H/SS_p_T*exp(p_H-p_T) )^(1-PP_Cons_ElasHF) ...
%             +(1-PP_Cons_FracH)*( SS_p_F/SS_p_T*exp(p_F-p_T) )^(1-PP_Cons_ElasHF);
e(id)   =   -p_T ...
            +PP_Cons_FracH*( p_H )+(1-PP_Cons_FracH)*( p_F );
% Inflation  (3)     
id      =   id+1;
e(id)   =   -(1+SS_Pir_N+Pir_N)+exp(p_N-Lag_p_N)*(1+SS_Pir+Pir);
id      =   id+1;
e(id)   =   -(1+SS_Pir_H+Pir_H)+exp(p_H-Lag_p_H)*(1+SS_Pir+Pir);
id      =   id+1;
e(id)   =   -(1+SS_Pir_F+Pir_F)+exp(p_F-Lag_p_F)*(1+SS_Pir+Pir); 
%--------------------------------------------------------------------------
% Government (5)
%-------------------------------------------------------------------------- 
% Monetary Policy (1)
id      =   id+1;
e(id)   =   -ir ...
            +PP_Taylor_ir*Lag_ir ...
            +PP_Taylor_Pi*Pir+PP_Taylor_Pi_N*Pir_N+PP_Taylor_Pi_H*Pir_H ...
            +PP_Taylor_dE*dE + PP_Taylor_dY*dY ...
            +Eta_M_dom;
% Fiscal Policy (4)
id      =   id+1;
e(id)   =   -( SS_T+T + (SS_B+Lag_B)/(1+SS_Pir+Pir) + SS_G+G ) ...
            +(SS_TAU+TAU)*( SS_w_N*SS_L_N*exp(w_N+L_N)+SS_w_H*SS_L_H*exp(w_H+L_H) ) ...
            +(SS_B+B)/(1+SS_ir+ir);
id      =   id+1;
e(id)   =   (1-PP_FiscalFlag_B) * ( -SS_B/(1+SS_Pir) + (SS_B+B)/(1+SS_Pir+PP_FiscalFlag_B_Real*Pirp) ) ...
            +(1-PP_FiscalFlag_G)*G;
id      =   id+1;
e(id)   =   (1-PP_FiscalFlag_T) * T + (1-PP_FiscalFlag_G)*G;
id      =   id+1;
e(id)   =   (1-PP_FiscalFlag_TAU) * TAU + (1-PP_FiscalFlag_G)*G;
%--------------------------------------------------------------------------
% Financial Market Returns (4)
%--------------------------------------------------------------------------
id      =   id+1;
e(id)   =   -ir_dom+ir;
id      =   id+1;
e(id)   =   -ir_ext+ir_s+(1-PP_Flag_dE)*dEp;
id      =   id+1;
e(id)   =   -rr_dom+ir_dom-Pirp;
id      =   id+1;
e(id)   =   -rr_ext+ir_s+dEp-Pirp;
%--------------------------------------------------------------------------
% Rest of the World (5)
%--------------------------------------------------------------------------
% Foreign Monetary Policy
id      =   id+1;
e(id)   =   -ir_s+Eta_M_ext;
% Demand for Tradeable Home Good
id      =   id+1;
e(id)   =   -SS_Y_H*exp(Y_H) ...
            +SS_C_H*exp(C_H) ...
            +exp( -PP_Cons_ElasHF*(p_H-p_F) )*( SS_C_H_s+Eta_Y_H*SS_Y_H );
% e(id)   =   -SS_Y_H*exp(Y_H) ...
%             +exp(Eta_Y_H)*SS_Y_H;
% e(id)   =   -Pir_H+Pir_H_s+dE;
% Law of One Price (1)
id      =   id+1;
e(id)   =   -Pir_F+Pir_F_s+dE;
% Foreign Inflation (2)
id      =   id+1;
e(id)   =   -Pir_H_s;
id      =   id+1;
e(id)   =   -Pir_F_s+Eta_p_F-Lag_Eta_p_F;
%--------------------------------------------------------------------------
% Market Clearing Conditions (4)
%--------------------------------------------------------------------------
id      =   id+1;
e(id)   =   -(SS_Bp_dom+Bp_dom)+SS_B+B;
id      =   id+1;
e(id)   =   -SS_Y_N*exp(Y_N)+SS_C_N*exp(C_N);
id      =   id+1;
e(id)   =   -SS_N_N*exp(N_N)+SS_L_N*exp(L_N);
id      =   id+1;
e(id)   =   -SS_N_H*exp(N_H)+SS_L_H*exp(L_H);

%--------------------------------------------------------------------------
% Accounting Equalities (3)
%--------------------------------------------------------------------------
id      =   id+1;
e(id)   =   -SS_Y*exp(Y) ...
            +SS_p_N*SS_Y_N*exp(p_N+Y_N)+SS_p_H*SS_Y_H*exp(p_H+Y_H);
id      =   id+1;
e(id)   =   -(SS_PAC+PAC) ...
            +(SS_PAC_N+PAC_N)+(SS_PAC_H+PAC_H);    
id      =   id+1;
e(id)   =   -(SS_Profit+Profit) ...
            +(SS_Profit_N+Profit_N)+(SS_Profit_H+Profit_H);

%--------------------------------------------------------------------------
% Lagged Variables (9)
%--------------------------------------------------------------------------
id      =   id+1;
e(id)   =   -Lag_p_Np+p_N;
id      =   id+1;
e(id)   =   -Lag_p_Tp+p_T;
id      =   id+1;
e(id)   =   -Lag_p_Hp+p_H;
id      =   id+1;
e(id)   =   -Lag_p_Fp+p_F;
id      =   id+1;
e(id)   =   -Lag_Bp+B;
id      =   id+1;
e(id)   =   -Lag_irp+ir;
id      =   id+1;
e(id)   =   -Lag_Eta_p_Fp+Eta_p_F;
id      =   id+1;
e(id)   =   -Lag_Yp+Y;

id      =   id+1;
e(id)   =   -dY + Y - Lag_Y;
%% Print the File
fid  =  fopen('Equ_NK.m','w+');
% Print the function head
fprintf(fid,'function [varargout]=Equ_NK(PP,SS,order)\n\n');

fprintf(fid,'%%');
for i=1:length(varlist)
    if mod(i,8)==0
        fprintf(fid,['''',char(varlist(i)),'''',',...\n %%']);
    else
        fprintf(fid,['''',char(varlist(i)),'''',',']);
    end
end

fprintf(fid,'\n \n ');
% Print the Structural Parameters
for i=1:length(pplist)
    sp  =   char(pplist(i));
    fprintf(fid,[sp,'=',strrep(sp,'PP_','PP.'),';\n']);
end
fprintf(fid,'\n ');
% Print the Steady State Values
for i=1:length(sslist)
    ss  =   char(sslist(i));
    fprintf(fid,[ss,'=',strrep(ss,'SS_','SS.'),';\n']);
end
fprintf(fid,'\n ');
% Print the Input Deviation Values
for i=1:length(varlist)
    fprintf(fid,[char(varlist(i)),'=0;\n']);
end
% Print the 0-order Equation
fprintf(fid,'if order==0\n');
fprintf(fid,['res=ones(',num2str(length(e)),',1);\n']);
for i=1:length(e)
    fprintf(fid,['res(',num2str(i),')=',char(e(i)),';\n']);
end
fprintf(fid,'varargout{1}=res;\n');
% Print the 1-order Equations
fprintf(fid,'elseif order==1\n');
for i=1:length(varlist)
    fprintf(fid,['fir_ord.',char(varlist(i)),'=[']);
    for j=1:length(e)
        der     =   jacobian(e(j),varlist(i));
        fprintf(fid,char(der));
        if j<length(e)
            fprintf(fid,';');
        end
    end
    fprintf(fid,'];\n');
end
fprintf(fid,'varargout{1}=fir_ord;\n');

fprintf(fid,'\n end\n');

fclose(fid);